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Abstract 


In order to study water-gas transport processes in the gas-diffusion-layer (GDL) of a proton exchange membrane (PEM) fuel cell system, a 
multiphase, multiple-relaxation-time lattice Boltzmann model is presented in this work. The model is based on the mean-field diffuse interface 
theory and can handle the multiphase flows with large density ratios and various viscosities. By using the standard bounce back boundary condition 
and an approximate average scheme for the non-slip and wetting boundary walls, respectively, detailed liquid-gas transportation in the GDL, 
in which exact boundary condition is difficult to be implemented, can be simulated. Unlike most of lattice Boltzmann methods based on the 
Bhatnagar—Gross—Krook collision operator, the present model shows a viscosity-independent velocity field, which is very important in simulating 
multiphase flows where various viscosities coexist. We validate our model by simulating a static droplet on a wetting wall and compare with 
theoretical predictions. Then, we simulate a water-gas flow in the GDL of a PEM fuel cell and investigate the saturation-dependent transport 


properties under different conditions. The results are shown to be qualitatively consistent with the previous numerical and theoretical works. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The proton exchange membrane (PEM) fuel cell is consid- 
ered as one of the best approaches in the search of new energy 
sources. In the PEM fuel cell, two reactant gases, hydrogen and 
oxygen, combine at a membrane of about 50 um thickness cov- 
ered with a catalytic layer. The membrane is surrounded by a 
gas-diffusion-layer (GDL) of about 200 wm thickness. In this 
setup, the GDL plays an important role because it has several 
specific functions such as providing a continuous transport of 
the reactant gases and electronic conductivity between anode 
and cathode [1]. However, since the PEM fuel cell is usually 
operated with humidified reactant gases at lower temperatures, 
a so-called “flooding” phenomenon, particularly on the cathode 
side where the water vapors condense and block the pores of 
the GDL, is possibly caused and hence degenerates the perfor- 
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mance of the fuel cell. Therefore, it is vital to investigate the 
liquid-gas transport in the GDL and study the transport proper- 
ties so that the diffusion media with optimal performance can 
be found. 

In the literature, to simulate liquid-gas two-phase flows in 
the GDL, Darcy’s law is usually applied for both phases with 
properties especially for two-phase effects. The properties are 
capillary pressure, which is the pressure difference between two 
phases, and relative permeability, which adjust the permeability 
of the Darcy’s law for each phase. These properties can be func- 
tions of liquid saturation, which is ratio of liquid volume to pore 
volume. Since the GDL is usually a complex random microstruc- 
ture, it is very difficult to characterize its exact geometry and 
so to calculate the transport properties [2]. Such difficulties 
can be illustrated in Fig. 1, which shows a micrograph of the 
Toray TGP-H-60 as a typical example and visualizations of the 
reconstructed three-dimensional (3D) geometries via the virtual 
material design and the dissipative particle dynamics (DPD), 
respectively. Obviously, one can clearly see the strong anisotropy 
of the medium structure. 
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Fig. 1. Micrograph of the Toray TGP-H-60 and visualizations of reconstructed 3D-geometries generated by the virtual material design (middle) and the dissipative 


particle dynamics (right), respectively. 


Conventionally, to elucidate the two-phase transport phenom- 
ena in the GDL, various numerical models have been developed. 
The simplest approach is a single-phase model, in which gas and 
liquid are considered as a single-fluid mixture and thus share the 
same velocity field. Also, the interfacial tension effect is com- 
pletely ignored. In this case, the total amount of water can be 
obtained by solving a single equation without distinguishing 
water vapor from liquid water. Once the total water concentra- 
tion field is obtained, one may allow for the water concentration 
going beyond the saturation level, essentially assuming super 
saturation in the gas phase [3-5]. The more rigorous approach 
to liquid water transportation in the GDL is a reformulation 
model in which the classical two-phase models is reformulated 
into a single equation, and the interfacial tension effect and GDL 
wettability, essential for successful fuel cell operation, are fully 
accounted for [6-10]. This model can efficiently produce most 
of the GDL characteristics, but neglect the local structure details. 
Another approach is in terms of volume-averaged or pore-level 
models [11—16]. These models assume local interfacial equilib- 
rium, namely, electrical, chemical, and thermal equilibrium at 
the pore level. Conditions of validity of local interfacial equi- 
librium were carefully defined. All of the above-cited models 
are, strictly speaking, macroscopic models, although theoreti- 
cal inconsistency may exist in some works, and by nature any 
interface boundary between two fluid phases is mesoscopic. 
Therefore, it is desirable to develop a mesoscopic model to study 
the two-fluid flows. 

Recently, as a mesoscopic model, the lattice Boltzmann (LB) 
method has become a promising tool to simulate the multi- 
phase flows [17-23]. The LB method originated from the lattice 
gas cellular automata (LGCA) [24,25]. Due to its kinetic ori- 
gin, the LB method poses some features that are different 
from the macroscopic models. These features include pro- 
gramming simplicity, intrinsic parallelism and straightforward 
resolution of complex boundaries and multiple fluid species. 
The later two features are quite appealing in simulating mul- 
tiphase flows in the GDL-like structures. In these cases, the 
fluid-solid interface is approximated by a zigzag approach so 
that the standard bounce back boundary condition can be directly 
applied, which reverse the momentum of fluid particles col- 


liding with a solid boundary by mimicking the particle—solid 
interaction for non-slip boundary condition. However, the major- 
ity of the LB methods used in modeling multiphase flows are 
based on the Bhatnagar—Gross—Krook (BGK) collision model 
[26,27] and they often encounter problems such as numeri- 
cal instabilities and viscosity-dependent velocity field [23,28]. 
Since various viscosities are usually presented simultaneously 
in the two-phase flows, the above limits become critical in 
accuracy of numerical simulations. To overcome these draw- 
backs, a multiple-relaxation-time (MRT) LB model is developed 
[22,23,28]. It has been demonstrated that the MRTLB model has 
better stability and accuracy in simulating the multiphase flows 
than the BGK counterparts [22,23]. However, most of the MRT 
multiphase models are based on the single-component potential 
model [18], which mimic two-phase flows by adjusting the inter- 
action potentials and can only be effective for two-phase flows 
with small density ratios (<100). A two-distribution-function 
MRTLB model was recently proposed and it overcomes most 
limits of the previous BGK and MRT model but requires an 
implicit treatment of the interface tension [29]. 

In this paper, a MRTLB model for multiphase flows with large 
density ratios and various viscosities is developed and applied 
to simulate water-gas transportations in the GDL of a PEM fuel 
cell. The model is based on the diffuse interface theory [30,31]. 
To handle the large density ratios and various viscosities, two dis- 
tribution functions are employed in this model, and one of which 
is used for modeling velocity field with total density of different 
phases, and the other is for tracking the interface by including 
effects of density differences between different phases. Unlike 
the previous MRT models, in which the interface tension and 
wetting boundary condition must be calculated in advance based 
on the static droplet on the solid walls, the present model explic- 
itly incorporates both factors. To catch the boundary effects, 
besides the standard bounce back condition used to mimic the 
non-slip boundary condition, an approximated average approach 
based on the Taylor series expansion is presented to model the 
wetting boundaries. 

The rest of this paper is organized as follows. In Section 2, 
a multiphase MRTLB model is presented. We first give a brief 
introduction of the diffuse interface theory, and then present the 
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MRTLB equations followed by numerical implementations of 
the non-slip and wetting boundary conditions. Section 3 includes 
two parts. In the first part this model is validated by simulating a 
static droplet on a wetting wall and the results are compared with 
theoretical solutions. In the second part, numerical simulations 
of the liquid-gas flows in the tentative model of GDL are carried 
out and the transport properties are investigated. We close with 
a conclusion in Section 4. 


2. Multiphase, multiple-relaxation-time lattice 
Boltzmann model 


2.1. Diffuse interface theory 


Our model is based on the diffuse interface theory [29,30]. 
Here, we consider a flow with two fluids or phases such as gas 
and liquid which have mass densities of pg and p; and viscosities 
of ng and nı, respectively. The whole flow system is characterized 
by the mass density p = p1 + pg, the viscosity n =n + Ng and the 
local order parameter gy denoting the density difference of two 
phases. The thermodynamic behavior of the system can then be 
described by the following free energy functional 


r= fav (v+ 5I¥0") + fasws (1) 


where V is a control volume and S is the surface area of the wet- 
ting wall. The first integral on the right-hand side is the standard 
Ginzburg-Landau expression with a bulk free energy density y 
defined by [31] 


y=AW-1), (2) 


where the parameters k and A determine the interface tension 
o = 4,/2kA/3, and interface width £ = ./2k/A [32]. With this 
form of the bulk free energy, the system will contribute to two 
equilibrium states, — 1 and +1, corresponding low and high den- 
sity fields, respectively. 

The second term in Eq. (1) accounts for specific wetting inter- 
actions at the solid—fluid boundaries due to the surface energy 
Ws. Functional minimization of F requires 


Ys = kon - V9), (3) 


which is evaluated at the solid wall, where n is the surface normal 
pointing to fluids [33]. This condition leads to a static contact 
angle @ at a flat wall in the absence of flow such that the wetting 
potential y has 


y = —kn- Vg = 2sign (5 — 0) ios (1 — cosh ) v2 k, 
(4) 
where £ = arccos(sin? 0). Given the free energy functional of Eq. 


(1), the dynamic behavior of the two-phase fluid is governed by 
the Navier-Stokes equations, 
dou 


y tV puu=—V-P+7V°u+ Fe, (5) 


and the Cahn-Hilliard equation, 
—+V-gu= MV°u, (6) 


with u being the macroscopic fluid velocity, Fp the body force 
and M is the mobility. The gradient of the pressure tensor P and 
the chemical potential jz are written as [30-32] 


V- P = Voc + ug) — uyg, (7) 
and 
f] 
u= É- Evry, (8) 
ag 


respectively, and cs is the sound speed given later. The bulk 
pressure of the system can be calculated as 


aw 1 
p=pc +9 -y-k |V- -IVg ). (9) 
dy 2 


2.2. MRTLB method 


In the framework of the MRT lattice Boltzmann method with 
N discrete velocities [28,29], Eqs. (5) and (6) can be solved by 
the following two equations: 


f(r + e;ôt, t + dt) = ffr, t) — O-' A (mp, t) — m‘(r, t)) 
+ dt (1- 50°'4;0) G(r, 1), (10) 


g(r + edt, t + ôt) 
= g(r, ft) — Q7' Ag(m,(r, i= mer, t)), (11) 


where i= 1,2, ...,N, fand g are the respective vectors of the den- 
sity and order-parameter distribution functions at lattice location 
r and time t. Ag(q=fg) is the diagonal relaxation matrix give by 


-> SaN=1)s (12) 


and Q is a Nx N matrix which linearly transforms the distri- 
bution functions f and g to the velocity moments mp and mg, 
respectively with 


my =Q-f, (13) 
m,=Q-g (14) 


The body-face symbols f, g, m and G denote the dimensional 
column vectors and have 


f(r, t):=(fo, fis ---, fv—1)', 


m(r, t):=(M fo, M f1, -- 


Aq = diag(Se0, Sal, -- 


g(r, 1):=(go, 81, -+> 8-1)", 
“Mpa: 

m,(F, f):=(190, Mg1, -.., MgN—1)", 
G(r, t):=(Go, G1,..., G-1)'. 


Here T denotes the transpose operator, and G represents the body 
forcing, of which components are 

(e; — u) 
z (Fp + uVo), (5) 


S 


G; = wi 
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where j4Vg comes from the pressure tensor (see Eq. (7)) and 
w; is the weight coefficients dependent on the discrete velocity 
model. For the three-dimensional 19-discrete-velocity (D3Q19) 
model with cs = 1/ /3(6r = ôt = ô = 1), the discrete velocities 
and weighting coefficients are 


(0,0,0), i=0O 
eq, = (+1, 0, 0), (0, +1, 0), (0,0, +1), i=1-—-6 : 
(+1, +1, 0), (1,0, +1), (0, 


i 0 

a. f= 

3 
wi = l i=1-—6 (16) 
i 18° t= , 

eg 

—, i=7-18 

36 


and the 19 equilibrium moments and their corresponding relax- 
ation coefficients in Ag are given in Table 1. 

It should be noted that all the nonlinear velocity terms 
can be omitted in the equilibria of {mei} because the fourth 
order isotropic lattice tensor is not required in recovering the 
Cahn-Hilliard equation (Eq. (6)) [32]. The relaxation parame- 
ters chosen in this work is followed the analysis of the MRT 
model [34] to reduce the viscosity-dependence velocity field. 
The transformation matrix Q for the D3Q19 model is given in 
[34]. 

The local density p, the order parameter g and the momentum 
j are calculated, respectively as 


p= oh a7) 


Table 1 


p=} si (18) 


ô 
j:=(jx, Jy, jz) = je; + —poG. 19 
f=Ga jys je) 2 fie + 5a (19) 


It should be indicated that the density p defined in our 
work can be thought as a nominal density, the real local 
density should be obtained in terms of the order parameter 
as P(P)= Pet (Y— Pg)(P1 — Pe M(G1— Pg), where pı and gy 
are the upper and lower limits of the order parameter and 
can be determined by the Maxwell’s equal-area rule. The 
reason we use the nominal density in our work is that the 
value of it changes little in the whole flow field and thus the 
incompressible limit of the LBM is satisfied. On the other 
hand, the density difference is reflected in the order parameter 
ọ and its related chemical potential u, which acts as a force 
and is added in the flow field. The magnitude of this force is 
controlled by the surface tension and the interface thickness 
only. Therefore, large density difference cases can be handled 
easily and stability is not critical in this model. Besides this, 
one can note that all the elements in Q are integers and the 
inverse of the matrix Q can be calculated by Q7! = Q™/(Q-Q") 
because Q is an orthogonal matrix, therefore, the program can 
be efficiently coded by hand to reduce time-consuming matrix 
operations and thus enhancing computational efficiency in 
running the program. The main steps in the implementation of 
the present method can be summarized as the following: 


e Step 1: Compute moments {mi} and {m<i} from the fluid 
density and velocity (initially setting m fi = mg and m fj = 


mp). 


Equilibrium moments and their corresponding relaxation coefficients of f and g (pọ is the initial value of the density, tf=3n/(pôt) + 0.5 and tg =M/I +0.5 and F is 


a free parameter controlled the mobility) 


mo =) sp=0 
e Lo jj 1 
m4 = —1l1p +19 +19 TE 
fl g ce D ng 
e o llj-j 
Mpa = 3 9 SN = Spl 
f2 2 2 po po“ 
ms, = jx,y,z sp5,7=0 
e 2., (2—syi) 
M £4.6,8 ~ 3 Je yz S f4,6,8 = 88s) 
ma — 3i Hoe 
f9 Po 
me _ 3k oN fence 
f10 — 20 ‘f10 = Sfl 
j} = j} 
mpm = >> sp =syi 
po i 
B-J 
mi, = -2 ~ Sfl2 = Sf 
f12 2p0 
mý = ah S3 =5fl 
s5 o 
ma = at Sfi4 = Sf 
po 
ms = Bh Sf15 = Sfl 
7 po 
mis- =0 Sflo-18 = Sf 


eq _ ee 
Ma =P Sg0=0 
Tu 
eq _ = 
m= —30¢ + 19-7 Sei = a 
s & 
Tu 
eq _ ree 
Moy = oe >a r Sg2 = Sg] 
S 
eq e. = 
M 93.5.7 = wo z 593,5,7 =O 
20. (2 — S¢1) 
e 5 dzz Sg4,6,8 = 8 Z 
846,8 300 , (8 — S¢1) 
eq _ wie 
Meg =0 S29 =Sgl 
eq _ See 
Mig =0 S210 = Sg1 
eq _ BERTS 
ma =0 Sell =Sg1 
eq 
M12 = S912 =Sel 
eq _ = 
M13 =0 Sg13 =Sgl 
eq _ = 
Moig = 0 Sg14 =Sgl 
eq _ a 
Mois =0 Sg15 = Sg] 
eq = 
me16-18 = 0 Sg16-18 =Sg1 
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e Step 2: Calculation the equilibria {mo} and {mei} in the 
moment space. 

e Step 3: Compute the post-collision density distribution func- 
tions by Eqs. (10) and (11). 

e Step 4: Advecting the distribution function {f;} and {g;}. 


2.3. Boundary conditions and numerical evaluation of the 
derivatives 


One of the advantages of the LB method over the conventional 
Navier-Stokes solvers is the implementation of the boundary 
condition on the complex geometries like the GDL structure 
shown in Fig. 1. Theoretically, this kind of geometry can only be 
approximated by the zigzag grids (see Fig. 2) due to the difficul- 
ties of experimental measuring or mathematical describing the 
exact boundary positions. In this situation, the standard bounce 
back boundary condition of the LB method provides an efficient 
and simple way to model the non-slip condition in terms of the 
distribution functions {f;}. As shown in Fig. 2, the unknown 
distribution functions (black dash arrow) coming from the solid 
part can be simply set to the known values of their corresponding 
ones with opposite directions (blue solid arrow) and the actual 
non-slip boundary can be realized at the one-half grid spacing 
beyond the last fluid node [23]. 

Since the wetting boundary condition (Eq. (4)) is an equi- 
librium condition, it is appropriate to impose it through the 
equilibria {mg} [35]. From Eq. (4), to set the equilibrium bound- 
ary condition we usually need to calculate the normal vector of 
the boundary so that the order parameter g and its second deriva- 


Fig. 2. 2D sketch of implementations of the boundary conditions of complex 
geometry (fluid points: blue solid circle; solid points: black solid circle; the 
nearest and next-nearest fluid-side points: inside the green dash elliptic circle). 


tives V9 on the boundary nodes are obtained. However, as we 
stressed earlier, it is difficult for the geometries like the GDL. 
To avoid this difficulty, in this work we use an approximation 
strategy to evaluate them. In Fig. 2, the boundary point P (red 
solid circle) is surrounded by some fluid points (blue solid circle) 
and some solid points (black solid circle). To obtain the values 
of y and V7¢ on the point P, Taylor series expansions of ¢ on 
the nearest and next-nearest neighboring fluid-side points of the 
point P (see Fig. 2) are carried out, and for each expansion point 
we have op = gf — 60,9 + O(67), which can be further approx- 
imated by pp = pfı — 50,y. Therefore, on average, the order 


= 
| 
| | |_| 


A = > 


— 


t=0 t= 100006 t= 400006 
(a) 0 = 60° 
= 5 T ~ T i i 
| | 
e © /e@ 
t=0 t= 100006 t = 400006 
(b) 0 = 120° 


Fig. 3. Time variations of the spreading and shrinking droplet shape at 0 =60° and 120°, respectively. The isosurface represents the order parameter g =0. 
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parameter g and its second derivatives V? on the point P can 
be approximated as: 


~ devGr — S9ngr) _ Zale + 181 O) 


PP N N (20) 
Vp ~ 2% (pf = PP _ dOnQP) 
N 
_ len -or +1818) an 
N 7. 


where ô is the unit grid space and gf; is the order parameters on 
the nearest and next-nearest neighboring fluid-side points of the 
point B;/=1,2,..., N with N representing the total number of the 
nearest and next-nearest fluid-side points of the point P; y can 
be calculated from Eq. (4) directly. The implementation of Eqs. 
(20) and (21) is straightforward by using the lattice stencils, and 
one can judge and count the fluid-side points by the streaming 
directions of the discrete velocities whether pointing to fluid 
field or not. 

In the LB framework, the first order and second order deriva- 
tives of g in the fluid field can also be conventionally calculated 
by using the lattice stencils so that the discretization errors are 
reduced. In this work, the following stencil scheme is used: 


V= 350, wig(r + est, t + Deij 

ôt 
= —24ọ(r, t) + 36% wio + e;ôt, t + ôt) 
7 6(6t)? , 


where j= 1, 2, 3 denotes the dimensions. 


(22) 


V9 (23) 


3. Results and discussion 
3.1. Model validation 


To validate the present model, a static droplet held on a wall is 
simulated, and the contact angle is tested. A hemispheric droplet 
with radius R placed on a wall is considered and it is brought 
to the equilibrium state at rest under the effects of the wet- 
ting potential y imposed on wall. The computational domain 


180 
150 
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90 


60 


Contact angle 0 


30 


0 
-0.04 


-0.02 0 
Wetting potential y 


0.02 0.04 


Fig. 4. Comparison of calculated contact angles with the theoretical predictions 
(Eq. (4)) (triangle, theoretical prediction; circle, simulation results). 


Fig. 5. A model of 3D gas-diffusion-layer with edges parallel to the main flow 
directing from left to right. Slice shows the velocity vectors at a mid-plane 
parallel to the main flow direction (6p = 1.7 x 1073, = 120°, n/ne = 100). 


~ 1.0E-03 

S ---A -- BGK (51X51X51) Era 
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aS A 

Š 2.5E-04 a 
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g 

v 

Š 

Q 0.0E+00 


0 0.1 0.2 0.3 0.4 0.5 0.6 
Kinetic viscosity n/p 


Fig. 6. Comparison of the calculated absolute permeability of the GDL as a 
function of the fluid kinetic viscosity. 


is divided into a 51 x 51 x 51 cubic lattice with 5=1 and its 
edges parallel to the x, y and z coordinate axes. The radius 
of the droplet R is set 105 or 155. The other parameters used 
are p1/Pg = 1000, m/ng = 100, pg = 1, ng =0.1, o =7.86 x 1072, 
&=4.56 and I =107?. The contact angle is changed in the range 


K,wp 
K,ywp 
K,wp by Eq.(29) n=4 
K,ywp by Eq.(29) n=4 


K,wp 
K,nwpP 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 
Saywp 


Fig. 7. Comparison of the relative permeabilities obtained by the present model 
and van Genuchten function (Eq. (29)). 
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of 45° < 0 < 145°. Periodic boundary conditions are applied in 
both x and y directions and non-slip/wetting boundary condi- 
tions are used at top and bottom walls in z direction. For this 
problem, the equilibrium radius Re and contact angle 6, of the 
droplet can be calculated by the following formula [36]: 


Cc 2 A]; ’ ( ) 
6 = arctan | ———— (25) 
c 2(R h) ’ 


Flow direction 


t = 500006t 


t = 1000006t 


t = 15000006t 


where A is the height of the droplet and b is the contact length 
of the droplet connecting to the wall. 

Fig. 3(a) and (b) shows the time variations of the droplet 
on the hydrophilic (@ = 60°) and hydrophobic (0 = 120°) walls, 
respectively. It is seen from Fig. 3(a), the droplet spreads as the 
time passes and finally reaches to its equilibrium shape with 
the calculated contact angle 6, =59°. As shown in Fig. 3(b), the 
connecting area between the droplet and the wall reduces as the 
time passes, and the droplet reaches its equilibrium shape with 

c = 122°. Fig. 4 shows the variation of the calculated contact 
angle with the wetting potential y. The present results are in good 
agreement with the theoretical curve obtained by Eq. (4). The 


Fig. 8. Snapshots of the liquid-gas transport processes in the model of GDL at 6p = 1.7 x 107? (left row) and 3.4 x 107? (right row), respectively. Flow direction is 


indicated by the arrow on the top of the figure. 
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agreement between the theory and our simulation indicates that 
the present method handles the two-phase flows with wetting 
boundaries correctly. 


3.2. Liquid-gas flows in GDL 


For the works presented in this section, we fit a 3D-geometry 
to a carbon-fiber structure of the GDL with a porosity of 80.3% 
and characteristic pore scale 10 um (see Fig. 1) by the DPD 
method [37]. The carbon fiber made in the GDL is hydrophobic, 
and therefore in the following the gas phase is denoted as the 
wetting phase (WP) while the liquid phase as non-wetting phase 
(NWP) for brevity. The computational domain is discretized by 
51x51x51 or 101 x 101 x 101 cubic lattice with the same 
coordinate arrangement as in the case of Section 3.1. Periodic 
boundary conditions are used on all sides of the domain for the 
velocity field. For the order parameter, it is assumed to be equal 
to 1 (liquid phase) at one side and —1 (gas phase) at the opposite 
side and the other four sides are treated as periodic bound- 
aries. Unless otherwise indicated, other parameters used in the 
simulations of this section are same as those used in Section 3.1. 


3.2.1. Absolute permeability and comparison between the 
BGK and MRT model 

The mathematical basis of immiscible two-phase flows in a 
porous medium is Darcy’s law, written as 
ū = — Kj) VP, (26) 
where U denotes the average flow velocity and K represents the 
absolute permeability tensor of the porous medium. The compo- 
nent kj; of the permeability tensors represents the permeability 
of the material in the direction parallel to the jth coordinate axis 
if the main flow direction is parallel to the ith coordinate axis. 
For convenience, here we assume that the edges of the sam- 
ple are parallel to the Cartesian coordinate axes. As a direct 
consequence, the off-diagonal elements of K are equal to zero. 
Furthermore, if the pressure gradient (i.e., the main flow direc- 
tion) is parallel to the ith coordinate axis, Darcy’s law reduces 
to the scalar equation and we have 

nLüi 

ép 

where id; is the average flow speed, L the length of the medium 
edge which is parallel to the ith coordinate axis, and dp stands 
for the pressure difference across the medium in the ith direc- 
tion. The present situation is visualized in Fig. 5, which shows 
streamlines of the velocity field at 5p = 1.7 x 1073, @= 120° and 
mi/ng = 100. 

The absolute permeability of the GDL with respect to differ- 
ent fluid kinetic viscosities is also calculated. In this simulation, 
a constant pressure drop 5p=1.7 x 107° is applied to the 
single-phase flow through the GDL on the 51 x 51 x 51 and 
101 x 101 x 101 grids. To show the advantage of the MRT 
model over its BGK counterparts, the simulation using the BGK 
model is carried out at the same time. Fig. 6 shows the calcu- 
lated absolute permeabilities at various viscosities obtained by 


kii = 


l (27) 


the MRT and BGK models. It is clearly seen that the permeabil- 
ities obtained by the MRT model remain essentially constant for 
both grid sets when the viscosity changes. However, the results 
of the BGK simulation change significantly and increase with 
increasing viscosity. Therefore, compared to the BGK model, 
the present MRT model is more suitable for simulating mul- 
tiphase flow systems, in which different viscosities coexist. In 
addition, by using the MRT model for the multiphase flow sys- 
tems, computational efficiency can be enhanced by using coarser 
grid set while desired accuracy is achieved. Based on the above 
observations, all the following simulations use the 51 x 51 x 51 
grid set. 


3.2.2. Relative permeability and effects of pressure drop, 
wettability and viscosity ratio 

The relative permeability K, is defined by an extension of 
Darcy’s law to 


K,(Sa)K _ 
n 


ū=-— P, (28) 
where Sa (ratio of phase volume over pore volume in GDL) 
reflects that the permeability depends on the saturation. Using 
Eq. (28) we can calculate K,(Sa) for both wetting and non- 
wetting phases. Theoretically, the relative permeabilities can be 
written in terms of the saturations of each phase as the well 
known van Genuchten function 


K,(Sa) = Sa!/3(1 — (1 — Saj 70/700” (29) 


where the values of n are in the range of 2 < n < 5. Fig. 7 shows 
the comparison of the calculated K,(Sa) by the present model 
and Eq. (29). Seen from this figure, the results obtained by the 
present model generally agree with the solutions of Eq. (29). 
With the increase of the NWP saturation, and the WP permeabil- 
ity decreases while the NWP permeability increases. In order to 
make a deep insight of the relative permeability, in the following 
we first investigate the dependency of the relative permeability 
on the pressure drop dp. In this case, we fix the contact angle 
@= 120° and viscosity ratio 7/ng = 100, and examine the effects 
of two values of pressure drop 5p = 1.7 x 1073 and 3.4 x 1073, 
respectively. 


Kwp (6p = 1.7e-3) 
k K,ywp (Op = 1.7e-3) 
0.8 4 x p K,wp (6p =3.4e-3) 


K,xwe (dp = 3.4e-3) 


0.6 


Kywt 
Kynwp 


0 0.1 0.2 0.3 0.4 0.5 0.6 
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Fig. 9. Relative permeability vs. NWP saturation for liquid-gas flow in GDL at 
different 6p with 6 = 120°. 
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Fig. 8 shows snapshots of the liquid-gas transport processes 
in the GDL at three different time steps at p= 1.7 x 107? and 
3.4 x 1073, respectively. For clearly displaying the liquid trans- 
portation, the GDL structure is not shown in this figure and the 
following Fig. 11. As illustrated in Fig. 8, larger dp drives the 
NWP to move faster than the small ôp because larger driving 
force more easily overcomes the capillary resistance due to the 


Flow direction 


interface tension. The influence of dp on the relative permeabil- 
ity K, for both NWP and WP is shown in Fig. 10. As shown in 
Fig. 9, when dp is increased, the WP relative permeability does 
not change much for the range of the NWP saturation. However, 
the NWP permeability is increased slightly at the intermedi- 
ate NWP saturation. These observations are in agreement with 
previous numerical and experimental works [22,36]. 


t = 500006t 


= 1000006t 


t = 1500006t 


Fig. 10. Snapshots of the liquid-gas transport processes in the model of GDL at contact angles 0 = 105° (left row) and 120° (right row), respectively. Flow direction 


is indicated by the arrow on the top of the figure. 
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Fig. 11. Relative permeability vs. NWP saturation for liquid-gas flow in GDL 
at different contact angles with 5p = 1.7 x 107? and m/ng = 100. 


Next, the effects of the wettability on the relative perme- 
abilities at 5p=1.7 x 107? and m/ng=100 are studied. The 
wettability of the wall is controlled by tuning the contact angle 0. 
Here, two contact angles 0 = 105° and 120° are examined. Fig. 10 
shows snapshots of the liquid-gas transport processes with the 
two contact angles in the GDL at three different time steps. As 
illustrated in this figure, at the smaller contact angle 0 = 105°, i.e. 
a stronger wettability of the GDL, the NWP transport faster than 
that at 0 = 120°. The reason behind this phenomenon is obvious: 
the stronger the wettability of the GDL has, the larger the driven 
force imposed on the liquid, and hence the easier the liquid 
movement overcomes the resistance of the capillary force. This 
mechanism also leads to a higher NWP relative permeability in 
the strong-wettability GDL than in the weak-wettability one, and 
this finding can be observed in Fig. 11, which shows K, for both 
NWP and WP as a function of the NWP saturation at contact 
angles 0 = 105° and 120°, respectively. Again, the observations 
from Figs. 10 and 11 are in agreement with previous numerical 
and experimental works [22,36]. 

The viscosity-ratio effects on the relative permeability at 
6p =1.7 x 1077 and 6 = 120° are further studied in Fig. 12, which 
shows the relative permeabilities of the NWP and WP as a func- 


K,we (m/n 100) 
K,we (ni /Ne= 100) 
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Fig. 12. Relative permeability vs. NWP saturation for liquid-gas flow in GDL 
at different viscosity ratios with ôp = 1.7 x 107? and 8 = 120°. 


tion of the NWP saturation for two viscosity ratios, n1/ng = 10 
and nı/ng = 100. Consistent with the findings in the previous 
works [22,36], we also observe that a decreased relative perme- 
ability of the NWP is caused by an increased viscosity ratio, 
especially at intermediate saturation. This phenomenon can be 
attributed to the so-called “lubricating” effect on the liquid flow 
due to a gas film, which is caused by that, gas flowing in rel- 
atively small pores is strongly coupled to the liquid flowing in 
the larger regions of pore structure, and makes liquid experience 
an apparent hydraulic slip [36]. The greater the viscosity ratio, 
the greater the hydraulic slip becomes. On the other hand, as 
shown in Fig. 12, the relative permeability of the WP slightly 
increases, which indicates that the lubricating effect cause an 
increase of the WP average velocity and hence the WP relative 
permeability. 


4. Conclusions 


In this paper, a MRTLB model is presented and applied to 
simulate water-gas transportations in the GDL of a PEM fuel 
cell. The model is based on the diffuse interface theory, and 
employs two distributions so that multiphase flows with large 
density ratios and various viscosities can be handled. The inter- 
face tension and wetting boundary condition in the present model 
are straightforwardly implemented. To numerically realize the 
boundary conditions for the complicated structure like GDL, 
besides the standard bounce back condition used for the non- 
slip condition, an approximated average scheme based on the 
extrapolation method is derived to mimic wetting boundaries. 
The numerical validation of the static droplet on a wetting wall 
shows good agreement with the theoretical predictions. Com- 
pared to the conventional macroscopic multiphase models like 
volume-averaged or pore-level models, which cannot track the 
interface itself, the present MRTLB model tracks the phase 
interface automatically and has advantages in programming sim- 
plicity, intrinsic parallelism and straightforward resolution of 
complex boundaries. Furthermore, the present model can easily 
handle the three-dimensional cases as shown in this paper. 

The water-gas transportation in a 3D modeled GDL structure 
is simulated and the transport properties including absolute and 
relative permeabilities are calculated. The effects of the pressure 
drop, wettability and viscosity ratio on the relative permeability 
are also investigated. The obtained results show the expected 
trend as in the previous numerical and experimental works. This 
investigation demonstrates that the present multiphase MRTLB 
model is potentially a viable method for the flows in fuel 
cells. 
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